Operational daily evapotranspiration mapping at field scale based on SSEBop model and spatiotemporal fusion of multi-source remote sensing data

Accurate understanding of daily evapotranspiration (ET) at field scale is of great significance for agricultural water resources management. The operational simplified surface energy balance (SSEBop) model has been applied to estimate field scale ET with Landsat satellite imagery. However, there is still uncertainty in the ET time reconstruction for cloudy days based on limited clear days’ Landsat ET fraction (ETf) computed by SSEBop. The Moderate Resolution Imaging Spectroradiometer (MODIS) remote sensing data can provide daily surface observation over clear-sky areas. This paper presented an enhanced gap-filling scheme for the SSEBop ET model, which improved the temporal resolution of Landsat ETf through the spatio-temporal fusion with SSEBop MODIS ETf on clear days and increased the time reconstruction accuracy of field-scale ET. The results were validated with the eddy covariance (EC) measurements over cropland in northwestern China. It indicated that the improved scheme performed better than the original SSEBop Landsat approach in daily ET estimation, with higher Nash-Sutcliffe efficiency (NSE, 0.75 vs. 0.70), lower root mean square error (RMSE, 0.95 mm·d-1 vs. 1.05 mm·d-1), and percent bias (PBias, 16.5% vs. 25.0%). This fusion method reduced the proportion of deviation (13.3% vs. 25.5%) in the total errors and made the random error the main proportion, which can be reduced over time and space in regional ET estimation. It also evidently improved the underestimation of crop ET by the SSEBop Landsat scheme during irrigation before sowing and could more accurately describe the synergistic changes of soil moisture and cropland ET. The proposed MODIS and Landsat ETf fusion can significantly improve the accuracy of SSEBop in estimating field-scale ET.


Introduction
Evapotranspiration (ET) is the loss of water through soil evaporation and plant transpiration, and it connects the global water and energy cycles. Understanding the cropland water use in agriculture water management requires detailed information about field-scale daily crop water found that the fusion of Landsat-MODIS vegetation index yielded higher estimation accuracy of ET than the fusion of reflectance using the Priestley-Taylor model. This is because more fusion times may amplify errors from remote sensing data and fusion algorithms. There is reduced error propagation because of the fusion of a single vegetation index band. This study aims to show how well a Landsat-MODIS data fusion framework can improve the SSEBop estimated actual daily ET at the field scale. This methodology initially used the MODIS and Landsat images to estimate actual ET under clear-sky conditions. Second, the ET fraction was selected as the critical parameter for data fusion by the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) [21]. Finally, the results were assessed by the automatic weather systems (AWS) and eddy covariance (EC) measurements over cropland in northwestern China and compared with the SSEBop estimates using only Landsat data. The analysis of the results was also listed in the discussion session.

Study area and data
The Heihe River Basin (HRB) is the second-largest inland river basin in northwest China's arid region, with the yearly precipitation ranging from 100 to 250 mm. The research area is located in an oasis-desert zone in the HRB's middle reaches, with a range of around 200 km×200 km. Maize is the principal crop of the irrigated fields, which grows from April to October and consumes a big part of the HRB's scarce water resources. From November to March, there is no crop in the fields. Fig 1 shows the land-use type in 2015 [22]; this region is mainly covered by croplands, grasslands, barren areas, etc. The grasslands area is close to the upper reaches of HRB.
The eddy covariance (EC) and automatic weather systems were installed through the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) projects [23,24]. The released EC products provide the half-hourly sensible heat flux (H) and latent heat flux (LE). The automatic weather system observations include the surface net radiation (Rn), the soil heat flux (G), soil temperature, soil moisture, wind velocity, precipitation, air temperature, humidity, etc. We used the Averaging Soil Thermocouple Probe (TCAV) method to correct the soil heat flow plate measured G [25]. We solved the energy non-closure and filled the LE gaps by the Bowen ratio correction and mean diurnal variation methods [11]. Further, we used the Flux Footprint Prediction (FFP) model and time series EC data to calculate the footprint climatology as the source area [2]. Fig 1 also shows the source area for the EC tower, and the ET pixels within this area were averaged for validation.
This research's coarse remote sensing data include MODIS land surface temperature, reflectance, and albedo products from the Terra satellite with a resolution from 500m to 1000m. We unified the resolution of MODIS products to 500m through bilinear interpolation. The Landsat 8 Collection 2 Level 2 Science Product (L2SP) provides the finer surface temperature and reflectance with a resolution of 30m. The Google Earth Engine acquired ERA5 daily averages reanalysis data was adopted as regional meteorological inputs. It includes maximum and minimum air temperature, mean air temperature, wind speed, dew point temperature, and surface pressure.

Model description
The SSEBop model suggests that available net radiation primarily drives the surface energy balance process. Differences in LST can quantify a decline in ET due to water stress and other factors. SSEBop estimates actual evapotranspiration (ET a ) by multiplying the ET fraction (ET f ) and reference ET (ET r ) as follows: where the k max is the coefficient that scales the grass reference ET (ET r ) into the level of a maximum ET experienced by an aerodynamically rougher crop, and a recommended value for k max is 1.2. The SSEBop developed a new modification to estimate the solar radiation by assuming the "average-sky" condition and the daily net radiation (Rn24) can be computed from ERA5 minimum and maximum air temperature. The advantage of this method is that it does not need the observation of sunshine hours. The ET r is calculated using the FAO-56 recommended PM equation with the ERA5 meteorological parameters and the Rn24. The ET f is computed by the pixel-based land surface temperature T s , the estimated temperature T c at the idealized "cold-wet" limit, and the predefined temperature difference dT between the"hot-dry" and"wet-cold" surfaces for each pixel as Eq (2): For a more detailed introduction to the SSEBop model, please refer to the original literature [9,12,16,26,27]. The SSEBop model uses 9 to 15 Landsat images every year to simulate annual daily evapotranspiration (ET) at the field scale. The daily ET values for dates between overpass images are derived by the daily ET and its nearest respective overpass ET f . The lack of high-resolution images will lead to errors in ET f simulation for many days. In this research, we used the MODIS and Landsat data to calculate the daily MODIS ET f and Landsat ET f on the overpass days and then applied the STARFM data fusion method to expand the time of Landsat like ET f for clear-sky images. Because the fusion process needs to consider the values of surrounding pixels, for partial clear sky images, we used the bilinear interpolation method to resample MODIS ET f . This operation enables the valid ET f pixel values to be as many as possible to ensure that the ET estimation error caused by data loss can be more avoided in the time reconstruction process. Fig 2 shows the overview of the enhanced SSEBop model methodology for estimating daily ET at the field scale. The results generated only by MODIS or Landsat are referred to as "SSEBop MODIS ET" and "SSEBop Landsat ET" after this, respectively. The results estimated by the fusion of MODIS and Landsat ET f are called "SSEBop MODIS-Landsat Fusioned ET" in the following.

Evaluation method
We chose the correlation coefficient(r), coefficient of determination (R 2 ), root mean square error (RMSE), Nash-Sutcliffe efficiency (NSE), and percent bias (PBias) statistics to evaluate the model results. The correlation coefficient (r) and coefficient of determination (R 2 ) describe the degree of collinearity between simulated and measured data. The RMSE indicates error in the units of the constituent of interest. The NSE is a normalized statistic that determines the degree of residual variance vs. recorded data variance. The NSE value indicates how closely the observed vs. simulated data graphic fits the 1:1 line [28]. NSE is computed as shown in Eq (3): where M is the modeled ET data point, O is the EC observed ET, and O mean is the mean of observed data for the evaluated constituent. Percent bias (PBias) measures the average tendency of the simulated data to be larger or smaller than their observed counterparts and is computed as: It is essential to identify the sources of bias and remove or reduce the bias primarily when the model is used in water balance or crop water-consuming monitoring. The separation of the error into bias and random components provides further helpful information for evaluating the ET models [26]. Obviously, for the same RMSE, a model with no bias is superior to a model with a bias since the random errors that dominate the RMSE will tend to compensate over time or space, thus conserving volumetric estimates. The mean bias error and random error were investigated separately for irrigated fields in the research area: where MBE (mm) is the mean bias error, n is the number of paired data points. MSE is the mean squared error (mm 2 ), RMSE (mm) is calculated as the square root of MSE. The calculation of the square of the random error (MSEe, mm 2 ) as the difference between MSE and (MBE) 2 can be made by rearranging the Eq (7). MSEe is the mean square error of the random error term "e". The MSEe shows the variability of the error itself from the average error. Since the square error terms are additive, the relative contribution of the random error and the square of the MBE can be expressed as a percentage of the MSE, which is the total square error between the modeled and observed values.  [29], respectively. Therefore, the estimated daily net radiation had relatively high accuracy. The R 2 was 0.73, the root mean square error (RMSE) was 3.21 MJm -2 d -1 , and the NSE was 0.61. The negative PBias indicated that the assumption of "averaged-sky" slightly overestimated the daily net radiation. The negative PBias of daily net radiation may cause higher dT and higher ET f in Eq (2). According to Senay et al. [27], use of observed net radiation will improve the results; however, to establish the hot/dry boundary limit, the average-sky condition suffices for operational applications and its advantage for computational simplicity. As shown by the grey dots in Fig 3, the overestimation of net radiation was mainly concentrated on cloudy days. The downward shortwave radiation in cloudy days is more affected by clouds, and the estimation method based on "average-sky" conditions may be insufficient. Since this research mainly evaluated the improvement of the time reconstruction method on the SSEBop model, we still adopted the original net radiation calculation method proposed in SSEBop.

Clear days' SSEBop Landsat ET
We collected Landsat 8 overpass images on clear days covering the study area. There were 7, 11, and 8 images from 2013 to 2015, respectively. The difference between the two results only lies in remote sensing data, which showed that it was more reasonable and feasible to use high-resolution Landsat data to estimate fieldscale ET. It was also found that both the SSEBop MODIS and Landsat underestimated the actual ET. There may be two sources of this systematic error. One is the selection of the k max parameter. The variation of k max parameter usually ranges from 1 to 1.3. Local optimization can effectively improve the results. On the one hand, it may be the potential ET error caused by the uncertainty of ERA5 meteorological data. However, the primary purpose of this paper is to analyze whether the fusion of MODIS and Landsat can make the SSEBop model better estimate the daily field-scale ET.

PLOS ONE
Estimating daily field scale evapotranspiration based on SSEBop and spatiotemporal data fusion different phenological periods of crops. From April to September, the ET of crops from sowing to harvest showed a trend of gradually increasing and then decreasing progressively, roughly within 2-8mm. In July and August, the daily ET reached 6-8mm in the peak season of crop growth. The precipitation in this area mainly occurred from June to August, and the daily rainfall changed within 0-15mm. The soil volumetric water content varied from 0.2 to 0.4 throughout the crop growing season.
It can be seen from the shadow of Fig 5 that there was a sudden peak of ET in March every year, and the precipitation was almost zero at this time. By further comparing the soil moisture, the soil moisture also showed the same increasing trend. Therefore, we can conclude that this was the irrigation time in spring before sowing. Although there was no crop growth in the field, irrigation in spring can suddenly increase the ET from nearly 0 to 4 mm/day and gradually decrease after sowing. The SSEBop MODIS ET and SSEBop MODIS-Landsat Fusioned ET can effectively reflect this change. At the same time, SSEBop Landsat ET is not ideal, which underestimated the increase of ET caused by irrigation. At this time, only considering NDVI is not enough. The surface temperature in the SSEBop model can effectively reflect the surface water status. The increase of soil moisture reduces the surface temperature, which increases ET f and the estimated ET. However, it requires the input of surface temperature at this critical time. The SSEBop Landsat only uses limited Landsat images on clear days, so it can not guarantee to obtain the change of land surface temperature information. It also proves that the method in this study can improve the simulation accuracy of the ET time scale. Table 1 shows the overall performance statistics of daily and monthly ET estimations using different methods. The NSE value close to 1 represents the overall high simulation accuracy of   25.0%) at daily scale and the monthly statistical results had a similar performance. Fig 6 shows the scatter plots for the three methods' ET estimation versus EC observations. Generally, the clear days' results had a higher coefficient of determination (R 2 ) than cloudy days. However, there was a lower RMSE on cloudy days because the reduced net radiation made the crop have higher canopy resistance, and the ET declined. From 2013 to 2015, the average daily ET on clear days observed by EC was 2.7 mm/day, while cloudy days was 1.7 mm/day. Through comparison, it was found that the SSEBop MODIS-Landsat Fusioned ET had better performance on clear and cloudy days than the SSEBop Landsat ET. Table 2 shows the error statics for the ET from the three methods. The relative contribution of the random error MSEe and the square of the MBE can be stated as a percentage of the MSE, which is the total square error between the modeled and observed values. The numbers in brackets in Table 2 represent the percentage of the corresponding error in the total error. From the statistical results, the MSE of SSEBop MODIS ET was the smallest, followed by SSE-Bop MODIS-Landsat Fusioned ET, and SSEBop Landsat ET was the largest. The percentage of random error in the total error of the three methods also had the same variation law, which was 90.6%, 86.7%, and 74.5%, respectively. A method with no bias (MBE 2 ) is better than a model with a bias for the same RMSE (MSE) since the random errors (MSEe) that dominate the RMSE (MSE) will tend to compensate over time or space. Therefore, we can conclude that the fusion of MODIS and Landsat ET f can effectively reduce the total error of field-scale daily ET estimation and increase the proportion of random error. In this way, random error reduction will be more evident than that of the original SSEBop Landsat ET in estimating long-time series ET, and the error accumulation phenomenon will be significantly improved.

ET spatial patterns
The validation and comparison of the previous chapters were carried out at the site scale, and the comparison of different methods in the region is equally important. Fig 7 shows the spatial distribution of annual evapotranspiration calculated by three methods. The same ET scale showed that the spatial distribution of SSEBop-MODIS and SSEBop MODIS-Landsat Fusioned ET was closer, significantly higher than that estimated by SSEBop-Landsat. The annual ET showed that the cropland ET was considerably higher than the bare soil area. The area with the highest ET value was concentrated in the range of the Heihe River, which indicated that the SSEBop model could estimate the ET of different underlying surfaces. Through ET f fusion, the SSEBop MODIS-Landsat Fusioned ET was closer to that of SSEBop MODIS in site validation and on the regional scale.
We resampled the SSEBop MODIS ET to the same size as Landsat ET and plotted the cropland pixels' ET density scatter maps of the three methods (Fig 8). Fig 8 shows that the correlation coefficient (r) between MODIS-Landsat Fusioned ET and SSEBop Landsat ET was the highest. The correlation between SSEBop MODIS and SSEBop MODIS-Landsat Fusioned ET was higher than that between SSEBop MODIS and SSEBop Landsat ET. The red part of the density scatters diagram is the concentrated distribution area of the central ET values. Through comparison, it showed that in the whole region, the farmland yearly ET of SSEBop MODIS was concentrated in 600-800 mm, while that of SSEBop MODIS-Landsat Fusioned ET was focused on 500-700 mm, and that of SSEBop Landsat was the lowest, concentrated in 400-600 mm. In terms of site scale validation, SSEBop MODIS can effectively reflect the subtle temporal changes of farmland ET due to its higher time resolution surface temperature. By comparing ET spatial distribution, we can conclude that SSEBop MODIS-Landsat Fusioned ET also had higher accuracy in spatial distribution because it was closer to the results of SSEBop MODIS.   SSEBop MODIS-Landsat Fusioned ET can provide daily field water consumption, which provides a basis for dekadal (10 days) cumulative ET mapping. Fig 9 shows the map of dekadal accumulated SSEBop MODIS-Landsat Fusioned ET in 2015. Under the same ET mapping scale, the dekadal ET map reflected ET's apparent change in crop growth. The cropland ET in this area was more significant from July to August, peaked in late July, and gradually decreased. It can be found from the Fig 9 that the cropland ET in late August was higher than that in mid-August. Comparing the precipitation data in Fig 5 indicates that the precipitation in late August just showed a small peak, and the sufficient water supply made the ET larger accordingly. It also confirms the effective response of surface water change to remote sensing surface temperature, which was reflected in the SSEBop ET results.

Discussion
The SSEBop model has been validated in various underlying surfaces, and the best performance was observed for croplands in recent researches [30]. Because of its simplicity and operability, it can be applied from field to continental scales. As a simplified surface energy balance At present, the daily field-scale evapotranspiration estimation using SSEBop and Landsat still faces uncertainty. The main reason is the lack of sufficient effective ET f to describe the change of surface vegetation and water conditions. Affected by the revisit period and weather factors, less than 10 Landsat images can be obtained in a year, which has significant limitations on the simulation of severe evapotranspiration changes in arid areas of China. Although the It effectively improves earth observation data's temporal and spatial resolution by MODIS and Landsat data fusion. However, this method has not been applied to the SSEBop model. The validation of this paper suggests that the simulation accuracy of field-scale daily evapotranspiration can be significantly improved by fusing MODIS and Landsat ET f computed by SSEBop under clear-sky conditions. In many previous studies, daily high-resolution evapotranspiration estimation is often carried out by fusing surface temperature and surface reflectance to improve model input parameters' temporal and spatial resolution. However, the spatio-temporal fusion algorithm has errors, and the errors caused by more fusion times will eventually be propagated into the output of the ET model. The spatio-temporal fusion needs to consider the contribution of effective pixels around the predicted pixels, and it is usually assumed that there is no obvious land use change between the predicted image and the images used for fusion on clear days. Therefore, we can choose the data completely free from cloud pollution for fusion as far as possible, and combine more available remote sensing images to avoid the error caused by the large time difference between the images used for fusion on clear days and predicted images. We also did the land surface temperature and vegetation fusion to test the SSEBop model, but worse results were found than the ET f fusion-based method proposed in this study (not shown).
Some other model uncertainties can be found in this study, such as the systematic underestimation of ET by the three methods. In the traditional crop coefficient approach for ET estimation in irrigated agricultural areas, the product of k max and ET f is equivalent to crop coefficient. In a recent study, k max was taken as the maximum crop coefficient value if it was available in the look-up table, which was calibrated by field data [30,31]. Therefore, k max is also a very sensitive parameter. For example, its 10% estimation error can lead to a 10% error of evapotranspiration. Therefore, such systematic underestimation errors may come from the reference evapotranspiration estimated by ERA5 reanalysis data and the k max value that has not been locally corrected by EC observations or soil water balance methods.
Field-scale daily evapotranspiration mapping provides more possibilities for precision irrigation management. In the future application research of the SSEBop model, we can pay more attention to the calibration of key parameters of different crop types. At the same time, we can compare and analyze the impact of different meteorological data products, such as ERA5, CMADS, and GLDAS data [32,33] on the performance of the model; or add more thermal infrared remote sensing data, such as sentinel-3A product [34], to improve data fusion accuracy and the performance of the SSEBop ET.

Conclusions
The ET f is a key inversion parameter in the SSEBop model, which contains comprehensive information on underlying surface temperature (moisture information) and vegetation index (crop growth). In this paper, the ET f values calculated by SSEBop with MODIS and Landsat were spatio-temporal fusioned to obtain the Landsat scale clear days' ET f , which was further used to estimate the daily ET at the field scale.
The validation at the site scale showed that the SSEBop MODIS had the best estimation accuracy, while SSEBop with only Landsat had the worst result. The estimation accuracy of daily high-resolution ET can be effectively improved by integrating MODIS and Landsat ET f . Firstly, the results from the fusion scheme's statistical parameters were better than the original scheme, with higher NSE values, lower RMSE, and PBias values. Secondly, the fusion scheme can significantly improve the underestimation of crop ET in the research area in March by the SSEBop Landsat method. The increase of ET in this period is caused by irrigation before sowing. It is difficult to capture the change of surface temperature at this time only using limited Landsat data. Finally, this fusion scheme can reduce the proportion of deviation and increase the proportion of random error in the total errors of SSEBop daily ET. The random error can be reduced with time and space, and the estimation accuracy in long-time series ET estimation will improve.
In future research, multi-source remote sensing data, such as the MODIS Aqua, ASTER, and Sentinel data, can be introduced to obtain more effective ET f values to improve spatiotemporal accuracy. The SSEBop model parameters can be calibrated in combination with ground data to further improve the estimation accuracy of the model.